library(Glimma)
library(limma)
library(GlimmaV2)
data(lymphomaRNAseq)
rnaseq <- lymphomaRNAseq
# 7 samples
# each instance is a gene
# attribute is number of counts per sample
# data.frame(rnaseq$counts)
# add lane
groups <- data.frame(genotype=rnaseq$samples$group,
lane= as.character(c(rep(4,5),3,3)),
miscCont=c(rep(4000,5),300,250),
miscDisc=c("blue","red",rep("green",5)))
# add libsize
groups <- rnaseq$samples$group
# fit
design <- model.matrix(~0+groups)
contrasts <- cbind(Smchd1null.vs.WT=c(-1,1))
# convert raw counts to logCPM values by automatically extracting libsizes and normalisation factors from x
vm <- voomWithQualityWeights(rnaseq, design=design)
fit <- lmFit(vm, design=design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)
dt <- decideTests(fit)
summary(dt)
## Smchd1null.vs.WT
## Down 1478
## NotSig 24424
## Up 1277
# log cpm is the attribute here
# averaged across all arrays in the original linear model fit
# each instance is a gene;
# the col vector names are the geneIDs
# the col vector attribute is mean expression
# data.frame(fit$Amean)
# for now assume fit object
# this is where info for the MD scatter plot comes from!
# data.frame(vm$E)
# fit$design
prepareXYData <- function(
x,
xvals=NULL,
yvals=NULL,
xlab=NULL,
ylab=NULL,
status=NULL)
{
# assume MA plot if no x/y are given
if (is.null(xvals) && is.null(yvals))
{
# for now assume fit object
xvals <- unname(x$Amean)
xlab <- "logcpm"
yvals <- unname(x$coefficients[,1])
ylab <- "logfc"
stopifnot(all(names(x$Amean) == names(x$coefficients[,1])))
table <- data.frame(xvals, yvals)
names(table) <- c(xlab, ylab)
table <- cbind(table, x$genes)
}
else
{
}
data <- list(x=xlab, y=ylab, table=table, cols=colnames(table))
return(list(plotType="XY", data=data))
}
rep(1,5)
## [1] 1 1 1 1 1
GlimmaV2("MA", fit, colour=dt, display.columns=c("PValue"))
## Warning in readLines(con): incomplete final line found on 'C:/Users/hkari/
## OneDrive/onenote/R/win-library/3.6/GlimmaV2/htmlwidgets/GlimmaV2.yaml'